#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include "mpzsmoothscan.h"

/*
    Copyright 2009-2010 Andrew V. Sutherland

    This file is part of "smoothrelation".

    smoothrelation is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 2 of the License, or
    (at your option) any later version.

    smoothrelation is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with smoothrelation.  If not, see <http://www.gnu.org/licenses/>.
*/

/*
    Given n integers a[i], computes o = \prod a[i], destroying a[] in the process.
    Copied from mpzutil.c
*/
void mpz_mul_set (mpz_t o, mpz_t a[], int n)
{
	register int i, j;
	
	if ( ! n ) { mpz_set_ui(o, 1); return; }
	while ( n >= 2 ) {
		for ( i = 0 ; i < (n+1)/2 ; i++ ) {
			if ( 2*i+1 < n ) {
				mpz_mul (a[i], a[2*i], a[2*i+1]);
			} else {
				mpz_set (a[i], a[2*i]);
			}
		}
		n = i;
	}
	if ( n == 2 ) {
		mpz_mul (o, a[0], a[1]);
	} else {
		mpz_set (o, a[0]);
	}
}


/*
	Given n integers a[i] each bounded by 2^{2^e}, finds the first
	integer in the list that is composed entirely of primes dividing P,
	or returns n if none are.  Uses Bernstein's algorithm.
	
	NOTE: The list size n must be a power of 2 and the parameter
	t should point to a preallocated list 2n entries.
        This function does not dynamically allocate any memory.

*/
int mpz_smoothscan (mpz_t *a, int n, int e, mpz_t P, mpz_t *t)
{
	int levels[64], lens[64];
	register int h, i, j, k, m, flag;
	levels[0] = 0;  lens[0] = n;
	for ( i = 1 ; lens[i-1] > 1 ; i++ ) { levels[i] = levels[i-1] + lens[i-1];  lens[i] = lens[i-1]>>1; }
	h = i-1;
	
	// copy list (this could be avoided but it is convenient to do so)
	for ( i = 0 ; i < n ; i++ ) mpz_set(t[i],a[i]);

	// build parital product tree of t[] (all nodes <= P)
	for ( k = 0 ; k < h ; k++ ) {
		flag = 0;
		for ( i = levels[k], j=levels[k+1] ; i < levels[k]+lens[k] ; i += 2, j++ ) {
			mpz_mul(t[j],t[i],t[i+1]);
			if ( mpz_cmp(t[j],P) > 0 ) { mpz_set(t[j],P); } else { flag = 1; }
		}
		if ( ! flag ) break;
	}
	if ( flag ) {
		if ( k != h ) { printf("weirdness in smoothscan, didn't reach the tree top\n"); exit (0); }
		mpz_mod(t[levels[h]],P,t[levels[h]]);
		k--;
	}

	// compute remainder tree
	for ( ; k >= 0 ; k-- ) {
		for ( i = levels[k], j = levels[k+1] ; i < levels[k]+lens[k] ; i+=2, j++ ) {
			// note, we could have both t child and parent equal to P, and we don't reduce in this case
			// if t child is equal to P, the parent must also be P, so we just test the child
			if ( mpz_cmp(t[i],P) ) mpz_mod(t[i],t[j],t[i]);
			if ( mpz_cmp(t[i+1],P) ) mpz_mod(t[i+1],t[j],t[i+1]);
		}
	}
	for ( i = 0 ; i < n ; i++ ) {
		if ( ! mpz_sgn(t[i]) ) return i;
		for ( j = 0 ; j < e ; j++ ) {
			mpz_mul(t[i],t[i],t[i]); 
			mpz_mod(t[i],t[i],a[i]);
			if ( ! mpz_sgn(t[i]) ) return i;
		}
	}
	return n;
}
